ARIMA 시계열 모형과 회귀모형을 결합할 경우 흥미로운 결과를 도출할 수 있다.
수식으로 표현하면 다음과 같다.
\[y_i = \hat{\beta} x_i + \epsilon_i \\ \epsilon_i \sim ARIMA(p,d,q)(P,D,Q)[m] \]
여기서 \(\epsilon_i\)를 계절 \(ARIMA(p,d,q)(P,D,Q)[m]\) 시계열 모형으로 식별한다.
UCI Machine Learning Repository: Air Quality Data Set를 상기 모형을 예측하는 데이터로 활용한다. 이탈리아 도시에서 IoT 센서를 통해 공기오염을 측정한 데이터다.
시간별로 다양한 센서를 통해 대기상태를 측정하는 것은 물론 해당 시간의 온도, 습도도 측정한다.
https://archive.ics.uci.edu/ml/datasets/Air+quality
date_time : Date (DD/MM/YYYY) Time (HH.MM.SS)co_gt : True hourly averaged concentration CO in mg/m^3 (reference analyzer)pt08_s1_co : PT08.S1 (tin oxide) hourly averaged sensor response (nominally CO targeted)nmhc_gt : True hourly averaged overall Non Metanic HydroCarbons concentration in microg/m^3 (reference analyzer)c6h6_gt : True hourly averaged Benzene concentration in microg/m^3 (reference analyzer)pt08_s2_nmhc : PT08.S2 (titania) hourly averaged sensor response (nominally NMHC targeted)n_ox_gt : True hourly averaged NOx concentration in ppb (reference analyzer)pt08_s3_n_ox : PT08.S3 (tungsten oxide) hourly averaged sensor response (nominally NOx targeted)no2_gt : True hourly averaged NO2 concentration in microg/m^3 (reference analyzer)pt08_s4_no2 : PT08.S4 (tungsten oxide) hourly averaged sensor response (nominally NO2 targeted)pt08_s5_o3 : PT08.S5 (indium oxide) hourly averaged sensor response (nominally O3 targeted)t : Temperature in °Crh : Relative Humidity (%)ah : AH Absolute Humiditylibrary(tidyverse)
# air_dat <- read_csv("https://gist.githubusercontent.com/sachinsdate/2f4fb513f3a5ca367ae1866e3f5b8613/raw/df4f6395f35f58bb59999e8ff11336f5bc5d76b9/air_quality_uci_mod.csv")
air_dat <- read_delim("data/AirQualityUCI.csv", delim = ";")
air_df <- air_dat %>%
# 변수명 확인
janitor::clean_names() %>%
# 결측값이 포함된 행과 칼럼 제거
filter(!is.na(date)) %>%
select(-x16, -x17) %>%
# 날짜와 시간을 결합하여 날짜시각형으로 자료변환
tidyr::unite("date_time", date:time) %>%
mutate(date_time = lubridate::dmy_hms(date_time)) %>%
# 소수점이 ,로 구분된 문자형 변수를 숫자형 변수로 변환
mutate(across(.cols = where(is.character), .fns = ~ str_replace(.x, ",", ".") %>% as.numeric)) %>%
# 시간순으로 정렬
arrange(date_time)
air_df %>%
slice(1:100) %>%
reactable::reactable()불러들인 날짜 시각에 대한 데이터 정합성 확인을 선제적으로 진행한다. 먼저, 시작시각과 종료시각을 파악하자. 2004년 3월부터 2005년 4월까지 대략 1년이 조금 넘는 기간 다양한 센서를 활용하여 공기에 포함된 화학물질을 측정한 데이터로 확인된다.
library(timetk)
air_df %>%
summarise(시작시각 = min(date_time),
종료시각 = max(date_time))# A tibble: 1 x 2
시작시각 종료시각
<dttm> <dttm>
1 2004-03-10 18:00:00 2005-04-04 14:00:00
불러들인 데이터의 행의 갯수를 9357 파악하고 pad_by_time() 함수로 교차확인하여 결측값이 없음을 확인한다.
library(timetk)
air_df %>%
pad_by_time(.date_var = date_time, .by = "h", .pad_value = NA) %>%
nrow()[1] 9357
시간 단위로 데이터가 너무 많아 시각화를 하는데 시간이 많이 걸리기 때문에 일별로 \(\frac{1}{24}\) 비율로 압축하여 일별 센서 측정값을 시각화한다.
air_df %>%
pivot_longer(-date_time, names_to = "센서", values_to = "측정값") %>%
group_by(센서) %>%
# 각 센서별 일별 측정값을 평균
summarise_by_time(.date_var = date_time,
.by = "day",
일평균 = mean(측정값, na.rm = TRUE)) %>%
ungroup() %>%
# 각 센서별 일별 시계열 시각화
plot_time_series(.date_var = date_time,
.value = 일평균,
.facet_var = 센서,
.facet_ncol = 4)본격적인 회귀분석 예측모형을 구축하기에 앞서 시각화를 통한 탐색적 데이터 분석을 통해 데이터에 대한 이해를 높여간다.
시각적으로 보면 특정시간에 측정값이 이상값으로 확인이 된다.
air_df %>%
select(date_time, pt08_s4_no2, t, ah) %>%
pivot_longer(-date_time, names_to = "센서", values_to = "측정값") %>%
# 각 센서별 일별 시계열 시각화
plot_time_series(.date_var = date_time,
.value = 측정값,
.facet_var = 센서,
.facet_ncol = 3)alpha 값을 조절하여 이상치에 해당되는 특이점을 식별한다.
air_df %>%
select(date_time, pt08_s4_no2, t, ah) %>%
pivot_longer(-date_time, names_to = "센서", values_to = "측정값") %>%
# 각 센서별 일별 시계열 시각화
plot_anomaly_diagnostics(.date_var = date_time,
.value = 측정값,
.facet_var = 센서,
.facet_ncol = 3,
.alpha = 0.04
)